function V_new = penaltyIteration( V, Vs, Mhat, tol, theta )
%PENALTYITERATION Summary of this function goes here
%  This function is Penalty American Constraint Iteration
%  Returns the option value for next time step
%  theta = 0: Fully Implicit (Default)
%  theta = 1/2: Crank-Nicolson
if nargin < 5
    theta = 0; end

Large = 1/tol;
N = length(V);
Pii = Large*(V < Vs);
P = sparse(1:N, 1:N, Pii);
I = speye(N);
quit = 0;
V_old = V;
while quit==0
    V_new = (I + (1-theta)*Mhat + P)\((I - theta*Mhat)*V + P*Vs);
    checkTop = abs(V_new - V_old);
    checkBot = max(1, abs(V_new));
    check = max(checkTop./checkBot);
    if check < tol
        quit = 1;
    end
    V_old = V_new;
    Pii = Large*(V_new < Vs);
    P = sparse(1:N, 1:N, Pii);
end

end

